using System;
using System.Collections.Generic;

namespace ICP_MS.Plugin.HPLC;

public class Calibration
{
	private int maxLocation = 0;

	private decimal[,] GussArray = new decimal[100, 100];

	private int dimension = 2;

	private decimal[] solutionArray = new decimal[100];

	public decimal[] CalculateLinearCurve(double[] dataX, double[] dataY)
	{
		int N = dataX.Length;
		decimal[] x = new decimal[N];
		decimal[] y = new decimal[N];
		decimal[] resultArray = new decimal[10];
		decimal sumy = default(decimal);
		decimal c;
		decimal d;
		decimal sumx;
		decimal sumxx;
		decimal sumxy;
		decimal b = (c = (d = (sumx = (sumxx = (sumxy = sumy)))));
		for (int i = 0; i < N; i++)
		{
			x[i] = (decimal)dataX[i];
			y[i] = (decimal)dataY[i];
		}
		for (int i = 0; i < x.Length; i++)
		{
			sumx += x[i];
			sumxx += x[i] * x[i];
			sumxy += x[i] * y[i];
			sumy += y[i];
		}
		GussArray[1, 1] = sumxx;
		GussArray[1, 2] = sumx;
		GussArray[1, 3] = sumxy;
		GussArray[2, 1] = sumx;
		GussArray[2, 2] = N;
		GussArray[2, 3] = sumy;
		dimension = 2;
		CalculateGaussEquation(dimension);
		resultArray[0] = solutionArray[1];
		resultArray[1] = solutionArray[2];
		return resultArray;
	}

	private decimal[] LinearRegression(decimal[] dataX, decimal[] dataY)
	{
		decimal sumx = default(decimal);
		decimal sumxx = default(decimal);
		decimal sumxy = default(decimal);
		decimal sumy = default(decimal);
		for (int i = 0; i < dataX.Length; i++)
		{
			sumx += dataX[i];
			sumxx += dataX[i] * dataX[i];
			sumxy += dataX[i] * dataY[i];
			sumy += dataY[i];
		}
		decimal a = ((decimal)dataX.Length * sumxy - sumx * sumy) / ((decimal)dataX.Length * sumxx - sumx * sumx);
		decimal b = (sumy - a * sumx) / (decimal)dataX.Length;
		return new decimal[2] { a, b };
	}

	public decimal[] CalculateLinearCurve(List<double> dataX, List<double> dataY)
	{
		int N = dataX.Count;
		decimal[] X = new decimal[N];
		decimal[] Y = new decimal[N];
		for (int i = 0; i < N; i++)
		{
			X[i] = (decimal)dataX[i];
			Y[i] = (decimal)dataY[i];
		}
		return LinearRegression(X, Y);
	}

	public decimal[] CalculateOtherCurve(double[] dataX, double[] dataY, CalibrationCurveType curveType)
	{
		int N = dataX.Length;
		decimal a = default(decimal);
		decimal b = default(decimal);
		decimal sumx = default(decimal);
		decimal sumxx = default(decimal);
		decimal sumxy = default(decimal);
		decimal sumy = default(decimal);
		decimal sumxxy = default(decimal);
		decimal sumxxx = default(decimal);
		decimal sumxxxx = default(decimal);
		decimal sumxxxxxx = default(decimal);
		decimal sumxxxxx = default(decimal);
		decimal sumxxxy = default(decimal);
		decimal[] x = new decimal[N];
		decimal[] y = new decimal[N];
		decimal[] p = new decimal[N];
		decimal[] q = new decimal[N];
		decimal[] resultArray = new decimal[10];
		for (int i = 0; i < N; i++)
		{
			x[i] = (decimal)dataX[i];
			y[i] = (decimal)dataY[i];
		}
		for (int i = 0; i < N; i++)
		{
			p[i] = (decimal)Math.Log((double)x[i]);
			q[i] = (decimal)Math.Log((double)y[i]);
		}
		decimal sumpp = default(decimal);
		decimal sump = default(decimal);
		decimal sumpq = default(decimal);
		decimal sumq = default(decimal);
		decimal sumpy = default(decimal);
		decimal sumxq = default(decimal);
		for (int i = 0; i < N; i++)
		{
			sump += p[i];
			sumpp += p[i] * p[i];
			sumpq += p[i] * q[i];
			sumq += q[i];
			sumpy += p[i] * y[i];
			sumxq += x[i] * q[i];
		}
		for (int i = 0; i < N; i++)
		{
			sumx += x[i];
			sumxx += x[i] * x[i];
			sumxxx += x[i] * x[i] * x[i];
			sumxxxx += x[i] * x[i] * x[i] * x[i];
			sumxy += x[i] * y[i];
			sumy += y[i];
			sumxxy += x[i] * x[i] * y[i];
			sumxxxxx += x[i] * x[i] * x[i] * x[i] * x[i];
			sumxxxxxx += x[i] * x[i] * x[i] * x[i] * x[i] * x[i];
			sumxxxy += x[i] * x[i] * x[i] * y[i];
		}
		switch (curveType)
		{
		case CalibrationCurveType.Average:
			solutionArray[1] = sumy / sumx;
			break;
		case CalibrationCurveType.Quadratic:
			GussArray[1, 1] = sumxxxx;
			GussArray[1, 2] = sumxxx;
			GussArray[1, 3] = sumxx;
			GussArray[1, 4] = sumxxy;
			GussArray[2, 1] = sumxxx;
			GussArray[2, 2] = sumxx;
			GussArray[2, 3] = sumx;
			GussArray[2, 4] = sumxy;
			GussArray[3, 1] = sumxx;
			GussArray[3, 2] = sumx;
			GussArray[3, 3] = N;
			GussArray[3, 4] = sumy;
			dimension = 3;
			CalculateGaussEquation(dimension);
			break;
		case CalibrationCurveType.Cubic:
			GussArray[1, 1] = sumxxxxxx;
			GussArray[1, 2] = sumxxxxx;
			GussArray[1, 3] = sumxxxx;
			GussArray[1, 4] = sumxxx;
			GussArray[1, 5] = sumxxxy;
			GussArray[2, 1] = sumxxxxx;
			GussArray[2, 2] = sumxxxx;
			GussArray[2, 3] = sumxxx;
			GussArray[2, 4] = sumxx;
			GussArray[2, 5] = sumxxy;
			GussArray[3, 1] = sumxxxx;
			GussArray[3, 2] = sumxxx;
			GussArray[3, 3] = sumxx;
			GussArray[3, 4] = sumx;
			GussArray[3, 5] = sumxy;
			GussArray[4, 1] = sumxxx;
			GussArray[4, 2] = sumxx;
			GussArray[4, 3] = sumx;
			GussArray[4, 4] = N;
			GussArray[4, 5] = sumy;
			dimension = 4;
			CalculateGaussEquation(dimension);
			break;
		case CalibrationCurveType.Exponential:
			GussArray[1, 1] = sumxx;
			GussArray[1, 2] = sumx;
			GussArray[1, 3] = sumxq;
			GussArray[2, 1] = sumx;
			GussArray[2, 2] = N;
			GussArray[2, 3] = sumq;
			dimension = 2;
			CalculateGaussEquation(dimension);
			a = (decimal)Math.Exp((double)solutionArray[2]);
			b = solutionArray[1];
			solutionArray[1] = a;
			solutionArray[2] = b;
			break;
		case CalibrationCurveType.Logarithm:
			GussArray[1, 1] = sumpp;
			GussArray[1, 2] = sump;
			GussArray[1, 3] = sumpy;
			GussArray[2, 1] = sump;
			GussArray[2, 2] = N;
			GussArray[2, 3] = sumy;
			dimension = 2;
			CalculateGaussEquation(dimension);
			a = solutionArray[1];
			b = (decimal)Math.Exp((double)solutionArray[2] / (double)solutionArray[1]);
			solutionArray[1] = a;
			solutionArray[2] = b;
			break;
		case CalibrationCurveType.Power:
			GussArray[1, 1] = sumpp;
			GussArray[1, 2] = sump;
			GussArray[1, 3] = sumpq;
			GussArray[2, 1] = sump;
			GussArray[2, 2] = N;
			GussArray[2, 3] = sumq;
			dimension = 2;
			CalculateGaussEquation(dimension);
			a = (decimal)Math.Exp((double)solutionArray[2]);
			b = solutionArray[1];
			solutionArray[1] = a;
			solutionArray[2] = b;
			break;
		}
		resultArray[0] = solutionArray[1];
		resultArray[1] = solutionArray[2];
		resultArray[2] = solutionArray[3];
		resultArray[3] = solutionArray[4];
		return resultArray;
	}

	public decimal[] CalculateOtherCurve(List<double> dataX, List<double> dataY, CalibrationCurveType curveType)
	{
		int N = dataX.Count;
		double[] X = new double[N];
		double[] Y = new double[N];
		dataX.CopyTo(X);
		dataY.CopyTo(Y);
		return CalculateOtherCurve(X, Y, curveType);
	}

	private decimal MaxInGuss(int k, int dim)
	{
		decimal temp = default(decimal);
		for (int i = k; i <= dim; i++)
		{
			if (Math.Abs(GussArray[i, k]) > temp)
			{
				temp = Math.Abs(GussArray[i, k]);
				maxLocation = i;
			}
		}
		return temp;
	}

	private void Exchange(int r, int k, int dim)
	{
		for (int i = 1; i <= dim + 1; i++)
		{
			GussArray[0, i] = GussArray[r, i];
		}
		for (int i = 1; i <= dim + 1; i++)
		{
			GussArray[r, i] = GussArray[k, i];
		}
		for (int i = 1; i <= dim + 1; i++)
		{
			GussArray[k, i] = GussArray[0, i];
		}
	}

	private void CalculateGaussEquation(int dim)
	{
		for (int m = 1; m <= dim - 1; m++)
		{
			decimal ark = MaxInGuss(m, dim);
			if (ark == 0m)
			{
				return;
			}
			if (maxLocation != m)
			{
				Exchange(maxLocation, m, dim);
			}
			for (int i = m + 1; i <= dim; i++)
			{
				for (int k = m + 1; k <= dim + 1; k++)
				{
					GussArray[i, k] -= GussArray[m, k] * GussArray[i, m] / GussArray[m, m];
				}
			}
		}
		solutionArray[dim] = GussArray[dim, dim + 1] / GussArray[dim, dim];
		for (int l = dim - 1; l >= 1; l--)
		{
			decimal me = default(decimal);
			for (int j = l + 1; j <= dim; j++)
			{
				me += GussArray[l, j] * solutionArray[j];
			}
			solutionArray[l] = (GussArray[l, dim + 1] - me) / GussArray[l, l];
		}
	}

	public decimal[] CalculateGaussCurve(List<double> dataX, List<double> dataY)
	{
		decimal a = default(decimal);
		decimal b = default(decimal);
		decimal c = default(decimal);
		List<double> y1 = new List<double>(dataY);
		List<double> x1 = new List<double>(dataX);
		for (int i = 0; i < y1.Count; i++)
		{
			y1[i] = Math.Log(y1[i]);
		}
		CalibrationCurveType cuveType = CalibrationCurveType.Quadratic;
		decimal[] abc = CalculateOtherCurve(x1.ToArray(), y1.ToArray(), cuveType);
		decimal a2 = abc[0];
		decimal b2 = abc[1];
		decimal c2 = abc[2];
		c = (decimal)Math.Sqrt(-1.0 / (double)a2);
		b = b2 * c * c / 2m;
		a = (decimal)Math.Exp((double)(c2 + b * b / c / c));
		abc[0] = a;
		abc[1] = b;
		abc[2] = c;
		return abc;
	}
}
